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Abstract. As a first step to meet the challenge to calculate the electronic structure 
and total energy of charged states of atoms and molecules adsorbed on ultrathin- 
insulating films supported by a metallic substrate using density functional theory 
(DFT), we have developed a simplified new DFT scheme that only describes the 
electrostatic interaction of an external charged system with a metal surface. This 
purely electrostatic interaction is obtained from the assumption that the electron 
densities of the two fragments (charged system and metal surface) are non-overlapping 
and neglecting non-local exchange correlation effects such as the van der Waals 
interactions between the two fragments. In addition, the response of the metal surface 
to the electrostatic potential from the charged system is treated to linear order, whereas 
the charged system is treated fully within DFT. In particular, we consider the classical 
perfect conductor (PC) model for the metal response, although our formalism is not 
limited to this approximation. The successful computational implementation of this 
new methodology and the PC model is exemplified by the case of a Na"^ cation outside 
a metal surface. 
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1. Introduction 

Recent progress in tlie studies of single atoms and molecules on ultrathin-insulating films 
supported by a metal substrate have opened up a new frontier in atomic scale science. 
Not only are these films able to decouple the electronic states of the adsorbates from 
the metal substrate states, but they also provide sufficient tunnelling current to be able 
to characterise, manipulate and image the adsorbed species by a scanning tunnelling 
microscope. Some important examples include charge state control of adsorbed species 
[H [2] , imaging of frontier orbitals |3l H] , coherent electron- nuclear coupling in molecular 
wires |5] and tunnelling- induced switching of adsorbed molecules P |7]. In particular, 
polar films such as NaCl bilayers are able to support multiple charge states of adsorbed 
species, which can be switched in a controlled manner by attachment and detachment of 
tunnelling electrons. Furthermore, these processes have been suggested to be responsible 
for the observed reversible bond formation in molecular switches [7j. To fully exploit 
the new potential opportunities provided by these systems, there is a need of theory 
to unravel the electronic and geometric structure and the excited state potentials of 
anionic and cationic states of the adsorbed species. 

Some important advances in this respect have been made by using Density 
Functional Theory (DFT) [H [9] calculations lU [HI lU 121 [7] . However, the delocalisation 
error in the current exchange-correlation functionals and the system size make the 
calculations very challenging. In fact, this error in current density functionals can 
result in fractional occupation of the adsorbed species and makes it difficult to identify 
various multiply charged states [10]. Even with these limitations, important progress 
has been made by using DFT+U [HI [12] to correct for the self-interaction behind the 
delocalisation error for a Ag atom [2] , despite it is generally a bit dubious for delocalised 
states such as for s states of atoms [13]. There has also been some developments in 
constraining the occupancy of adsorbate orbitals, but this requires that one has to 
identify the appropriate Kohn-Sham orbital, which can still be challenging due to the 
mixing of adsorbate and substrate states albeit being very small [H]. However, when 
considering, for example, organic molecules adsorbed on an ultrathin-insulating film 
supported by a metal surface, the number of metal atoms and the associated number of 
metal electrons makes the calculations to become prohibitively large. 

As a first step to try to overcome these limitations, we present a new simplified 
DFT method for the calculation of the total energy, ionic (Hellman-Feynman) forces 
and electronic structure of a charged system placed in front of a metal surface. In this 
method, the charged system is fully treated using DFT, the interaction between the two 
systems is assumed to be electrostatic and the density response of the metal surface to 
the electrostatic interaction is treated to linear order. The proposed method has two 
main advantages. Firstly, the computational effort is significantly decreased, since the 
metal electron states are not treated explicitly but only implicitly through their density 
response to an external electrostatic field, which can be captured using a simple model. 
Here, we have explicitly considered the classical perfect conductor (PC) model for the 
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metal response in which the screening length by the conduction electrons is assumed to 
be zero and the screening charge only resides on the perfect conductor plane. Secondly, 
different charge states of the charged system can be handled directly using this method. 

We have implemented this new methodology in the VASP code |23] using the simple 
PC model for the metal response. A critical test of this implementation is provided by 
the example of the Na"*" ion outside a perfect conductor. Here, the six electrons in the p 
semi-core of the Na"*" ion were included in the calculations. The results show an excellent 
agreement between the calculated Hellman-Feynman force on the ion and the negative 
gradient of the interaction energy even in the region close to the perfect conductor plane, 
where the electron density of the ion overlaps with the induced charge density at perfect 
conductor plane. Furthermore, the small polarisability of the ion makes the calculated 
interaction energy to be in close agreement with the analytical result for a point charge 
model outside the perfect conductor. 

In applying this method to study charged adsorbates deposited over ultrathin- 
insulating films supported by a metal substrate, we need to extend the PC model 
by including the interactions between the film and the metal substrate arising from 
overlapping electron densities and van der Waals interactions. We plan to capture these 
interactions through a simple parametrised force field between the metal substrate and 
the ions in the film, where parameters will be obtained by fitting the force field to DFT 
calculations of the forces between the film (without adsorbates) and the metal substrate. 
A successful implementation of this force field should enable us to describe, for example, 
the observed reversible bond formation of a molecular switch induced by a tunnelling 
electron and hole attachment to form ionic adsorbate states [7]. The presentation of 
this scheme and associated results are deferred to a future publication. 

This paper is organised as follows. In the next Section [21 we present the theoretical 
background of the proposed method. First, the effective energy density functional and 
the associated Hellman-Feynman forces for the charged system are derived in Subsection 
12.11 under the assumption that the electrons of the charged system and the metal 
surface are distinguishable. The total energy of the metal surface is expanded up to 
second order in the resulting electrostatic interaction between the two systems. In the 
next Subsection \2.2\ the non-trivial effects on the electrostatics when imposing periodic 
boundary conditions are analysed and appropriate dipole corrections for the total energy 
and the Hellman-Feynman forces are derived. The classical perfect conductor model for 
the metal surface response is then introduced in the next Subsection 12.31 Here, we also 
present analytical expressions for the interaction energy and force in the case of a point 
charge (Subsection 12.3. ip in the presence of periodic boundary conditions, and discuss 
possible refinements of the perfect conductor model (Subsection 12. 3. 2p . In Section [3l the 
keys steps in implementing the PC model in the plane wave code VASP are presented. As 
a first test of this methodology, we show the results of the computational implementation 
of the PC model for the case of Na"*" ion outside a PC (Section H]) and finally give some 
concluding remarks in Section [5l 
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2. Theoretical framework 

In this Section we develop a density functional theoretical description of a charged 
system placed in front of a metal surface based on the assumption that the electrons of 
both fragments are distinguishable, so that the interaction between the two systems is 
only electrostatic. The corresponding approximation becomes valid when the electron 
densities of the two systems are well separated. We analyse how the system can be 
treated in a supercell geometry and derive expressions for the dipole correction of the 
Kohn-Sham potential and the Hellmann-Feynman forces. As an explicit model for the 
metal surface, we focus on the classical perfect conductor (PC) model of the metal 
surface response, but also discuss possible extensions of this approximation. For the 
particular case of a point charge, we present explicit expressions for the interaction 
energy with a perfect conductor. Throughout the paper, we will make use of electrostatic 
units. 

2.1. Approximate energy density functional 

We consider a closed, charged system S of electrons and ions outside a metal surface M. 
In contrast to the system S, the metal surface will be treated as an open system, kept at 
a constant chemical potential fi. The electron densities of S and M are assumed to be 
non-overlapping and are given by ns{r) and ^^(r), respectively. In this case the kinetic 
energy functional of the non-interacting electrons is additive in the electron densities of 
S and M. Furthermore, we will neglect any non-local contributions between S and M 
such as the van der Waals interaction so that the exchange- correlation energy functional 
is also additive in these electron densities. The total energy density functional for the 
combined system is then given by the following expression: 



Here -Esfrig] and Em[nm] are the energy density functional for the isolated 5* and M, 
respectively, whereas pm(r) is the charge density of M and Nm the total number of 
electrons of M. The potential 0s (r) is defined as the electrostatic potential generated 
from the total charge density Ps(r) = — ens(r) + Pi{r) of electrons and ions in S and is 
given by. 



Note that the approximation in Eqn.([T]) is still meaningful for overlapping electron 
densities. Furthermore, any neglected non-additive contributions to the kinetic energy 
and the exchange-correlation potential including any non-local contribution between M 
and S, such as the van der Waals interaction, will be accounted by introducing a force 
field between M and S. 

An effective total energy density functional for S outside M in term of the electron 
density ^^(r) is obtained by minimising the energy density functional in Eqn.(IT]) with 




(1) 




(2) 
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respect to for fixed and /x, 

/i + e0s(r) (3) 



and, since /i is constant, we conclude that the ground state electron density of M in 
the presence of S* is a functional of n^, i. e. Um = nm[ns]- The change in energy of M 
induced by the presence of S is defined as, 

/\Err\ns] = {Em[nm] " Em[nmo\) " /i(A^m " A^mo) + J Pm(r)0s(r)rfr (4) 

where n^o and A^^o are the unperturbed ground state electronic density and total 
number of electrons of M, respectively. Expanding AEm[ns] up to second order in 
0s (r), one obtains 

The two terms in Eqn. ([5]) can now be cast in a more familiar form using the result 
( ^ ^^^^ ff^^ 

which follows from the fact that AEm[ns\-, defined in Eqn.dl]), is stationary with respect 
to variations in nm(r) for fixed 0s . Thus, the first term of the r.h.s in Eqn. ([5]) is 
determined by the unperturbed charge density p^o of the metal M, whereas the second 
term is given by the induced charge density Pm(i(r) in M, which, to linear order in 0^, 
is given by 

Aw(r) = j ^-^fs{r')dr' . (7) 

Note that the symmetry of the kernel in Eqn. ([5]) and Eqn. ([6]) give rise to the following 
important symmetry condition for the density response kernel in Eqn. ([7]), 

50s(rO 50s(r) ' ^ ' 

Finally, we obtain the following effective total energy functional for 5* outside M from 
Eqns.dS]) and ([7]) as, 

Eeff[ns] = Es[ns] + AEm[ns] = 

Es[ns] + J pmo(r)0s(r)(ir + i J pi„d(r)0s(r)rfr . (9) 

The interaction energy between S and M is then simply defined as 

Eint = E^fflus] - Es[nso] = (10) 
Es[ns] - Es[nso] + J p,„o(r)0s(r)c?r + ^ J Pi„rf(r)0s(r)c?r , (11) 

where n^o is the unperturbed ground state electron density of S. The contribution from 
the first two terms in Eqn. fllip gives the interaction energy from the polarisation of the 
system S, whereas the contributions from the third and fourth terms give the interaction 
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energies of S with the unperturbed surface density and the polarisation of the metal 
surface, respectively. 

The condition that the chemical potential of M should be constant imposes an 
important constraint on the electrostatic potential from the induced charge densities in 
M and S. From the functional derivative of Eqn.(l3]) it follows 

= T TT + T TT - e(l)m " e0,(r) , (12) 

dnm{r) dnm{r) 
where 0m (r) is the electrostatic potential 

0m(r) = (t>mo{r) + (pindir) (13) 

obtained from Pm(i") = Pmo(r) + Pind{^)- Inserting Eqn. (fT3!) in Eqn. (fT2!) . one finally has 



We now assume that the metal surface M is perpendicular to the z-axis and it is a 
semi- infinite system that extends to — cxd. The contribution from the first three terms of 
the r.h.s. of Eqn. ffT^ gives the chemical potential in the absence of the perturbation by 
5" and, since the electron density is unperturbed far inside the metal surface {z — > —00), 
this contribution is constant. Therefore, the contribution from the last two terms has 
to be zero in this region, 

0s(r) + (pind{r) -> , 2; -> -00 . (15) 

The next step is to determine how the Kohn-Sham (K-S) potential of S changes in 
the presence of M. Since the K-S potential is determined by the functional derivative 
of the electrostatic energy and exchange correlation energy terms with respect to the 
electron density, we need to investigate how these terms are modified. The exchange 
correlation energy term is a universal functional of the density and does not change. On 
the other hand, according to Eqn.([9]), the electrostatic energy term in Es[ns], is defined 
here as, 

^ I p.(r)0,(r)rfr , (16) 

and changes in the presence of M to the effective electrostatic energy 

E'J^iris] = E,i[ns] + J Pmo(r)0s(r)c/r + ^ J p,w(r)0.(r)cir . (17) 

Thus, only the electrostatic potential energy term — e0s(r) in the K-S potential is 
modified and is now given by, 

- e0:rK] = = -e[Mr) + <Pmo{r) + 0.„.(r)] . (18) 

In deriving this result from Eqns.(|5]), (E]) and ([7]) , we have made use of the symmetry 
condition in Eqn.(|8]). Therefore, the only change of the K-S potential in the presence of 
M is simply the inclusion of the electrostatic potential from the metal surface. 
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The presence of the metal surface will now change the forces F/ on the nuclei / 
at positions Rj and charges eZj in S, and this modification will be only through the 
electrostatic potential from the metal surface. Since the effective total energy functional 
is stationary around the ground state electron density and only the electrostatic term 
has an explicit dependence on the nuclei positions, the forces F/ are given by 

r A F"-^-^ 

?eff _ / '^^el 



Fj = -VjE:1^ = -J j^VnM^) dr. (19) 

In a similar manner as for the functional derivative with respect to the electron density 
in Eqn.( !T8|) . one obtains that 

-^ = 0,(r) + 0„(r). (20) 

The expression for the forces in Eqn. f|T9|) can now be cast in a more familiar Hellman- 
Feynman form using Eqn.( l20l) and, since 

VR,p,(r) = -eZ,V5(r-R,), (21) 

one obtains, after partial integration, the total force acting on the ion / 

F, = -eZ,V[0. + 0„](R/). (22) 

Therefore, the force on a nuclei is still given by the total electrostatic field at the nuclei, 
but now includes the electrostatic contribution from the metal surface. 

2.2. Periodic boundary conditions and dipole corrections 

Since the combined system will be represented in a finite supercell, we need to discuss 
the effects of imposing periodic boundary conditions (PBCs) [T5l [16] . Firstly, we 
will introduce PBCs in the lateral directions along the planar surface of the semi- 
infinite metal and state some general key properties of the behaviour of the electrostatic 
potentials. In addition, the system will have a net dipole moment component along the 
direction perpendicular to the metal surface and, when imposing PBCs in this direction, 
one has to introduce appropriate dipole corrections for the total energy, K-S potential 
and the forces on the nuclei as detailed in the following. 

In discussing the behaviour of the electrostatic fields and charge densities in the presence 
of lateral PBCs, it is convenient to introduce a plane wave representation over two- 
dimensional (2D) reciprocal lattice vectors G of the two dimensional surface unit cell A 
with area A. The 2D plane wave expansion of a charge density p(r) is here defined as, 

p(r) = 5^p(^,G)exp[zG.R] (23) 

G 

where the plane wave coefficients p{z, G) are given by 

p{z,G) = ^ p{K,z)exp[-iG.R]dR (24) 
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and (R, z) = r. The electrostatic potential from the charge density can now be obtained 
using the associated plane wave representation of the Coulomb kernel, 
1 271. „ 27r 



- + ^ exp[-G|;^-z'|]exp[zG.(R-R')] (25) 

^ I G^O 

where |G| = G. In particular, the laterally averaged potential, (p^z) = (p{z,G = 0) is 
determined by the laterally averaged density, p{z) = p{z, G = 0), as, 

0(2) = -27r J \z- z'\p{z)dz' + 00 • (26) 

Note that (t){z) is only determined up to a constant 0O; which has to be fixed by the 
boundary conditions as discussed later. Finally, the remaining laterally varying part of 
the potential, which has contributions only from the non-zero G plane wave coefficients, 
is given by 



J]^y" exp[-G'|2-/|]exp[iG.R]p(/,G)ci427) 



G^O 

and decays exponentially away from a localised charge distribution. 

In the case of a semi-infinite metal surface, the external electric field from the 
charged system sets up an electron current towards the surface and builds up a localised 
surface charge distribution and an electric field that opposes the external field. In a 
stationary situation, the perpendicular components of the electric field and current far 
inside the surface have to be zero. This electric field is determined by the laterally 
averaged electrostatic potential (pi^z — > —00) and, according to Eqn. (!26l) . it is given by 

E,{z) = -^{z) ^ 27r y p{z')dz' , z ^ -oc (28) 

Therefore, in order for the electric field to be zero far inside in the metal, the total 
charge of S and M has to be zero. Since the surface charge of the unperturbed metal is 
zero 

J^J p^o{r)dRdz = , (29) 

the total induced surface charge has to be equal in magnitude to the total charge Qs of 
S but with opposite sign, 

^ J Pind{r)dRdz = -Q, . (30) 

The energy of the vacuum level for the unperturbed metal is defined to be equal to zero, 
(0mo(^) — 0, 2 — i- 00), then 

(Pmoiz) = -271 J [\Z - Z'\ + z']pmo{z')dz' . (31) 

The undetermined constant 0o in 4>s{z) + (t>ind{z) can now be determined from the 
condition in Eqn.f llSp that the chemical potential of the metal is constant: 

4>s{z) + (pindiz) -^0,2;-^ -00 . (32) 
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Using Eqns.( 12B]) and (15^ . this condition gives, 

Mz) + ^rUz) = -2n [ [\z - z'\ - z'][ps{z') + pind{z')]dz' (33) 



Since the screening length in a metal is typically short on the order of 1 A, the 
surface charge distribution of the semi-infinite metal surface M is highly localised and 
the combined system is neutral, so that the total system can be confined in a supercell 
V of volume V by introducing a PBC in the perpendicular direction. However, the 
separation of charges between M and S gives rise to a dipole moment and a long-ranged 
potential (t){z) that has to be treated carefully. In addition, one also has to ensure 
that the length L of the supercell along the z-direction is large enough so that the 
short-ranged part of the potential 0'(r), Eqn.( [27l) . is confined within V . 

In order to proceed, we first need to discuss the electrostatics in a supercell. The 
electrostatic potential 0(r) from a charge distribution p(r) is most efficiently obtained 
by introducing a three dimensional (3D) plane wave representation of the densities and 
electrostatic potential over the 3D reciprocal lattice vectors g of the supercell. The 
plane wave expansion and coefficients p(g) of a density p(r) are defined here as, 

= EgP(g)exp[ig.r] ^^^^ 
= ^/vPWexp[-ig.r]rfr . 

The plane wave coefficients of the electrostatic field 0(g) from the density p(r) are 
obtained from the solution of the Poisson equation in reciprocal space and are given by. 




, g = 



9^ 

Note that this solution for (/)(r) from a p(r) with a net charge corresponds to the 
electrostatic potential from p(r) compensated with a uniform neutralising background, 
— p(g = 0) inside the supercell. Furthermore, the undetermined constant of 0(r) is set 
so that the average value 0(g = 0) of 0(r) is zero. 

In the case when the system has a net dipole moment, an artificial uniform electrical 
field has to be generated in the supercell in order to fulfil the PBCs. An efficient solution 
to this problem was first proposed by Neugebauer and Scheffler [17] and later corrected 
by Bengtsson [IB]. In this approach, a net dipole moment of a charge distribution p(r) 
(perpendicular to the z-axis) within the supercell is compensated by introducing a dipole 
layer aX z = z^ip, 

Pdip(r) = mS'{z - Zdip) (36) 
where the surface dipole moment density is defined as, 

fn = — [ p{r)zdr (37) 



and the resulting dipole potential in the supercell is given by 

'z 1 
L~2 



4>dip{z) = ATxm 



0<z <L . (38) 
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when Zdip = L. By correcting the electrostatic potential in the supercell using this dipole 
potential, one obtains the dipole-corrected potential 

0dipcor(r) = 0(r) + (pdipiz) . (39) 

This corrected potential in the supercell is now equal, up to a constant, to the 
electrostatic potential from S and M in the absence of the corresponding PBC in the 
direction perpendicular to the surface. 

The dipole-corrected effective electrostatic energy is obtained by correcting the 
electrostatic potentials 0s (r) + 0md(i") and 0mo(r) in E^f^ (Eqn. [TTj) as 

E'^Jf4^P ^ f p^(r)[0„o(r) + <f),,^o{z) + + (40) 
Jv 

^ J ps(r)[0s(r) + (pindir) + (l)dipi{z) + 4>i]dr, (41) 

where (pdipi{z) and (pdipoiz) are the dipole potentials from the surface dipole moment 
densities of Ps(r) + Pind(j) and Pm,o(r), respectively. The constants 0i and 0o have to be 
determined such that the conditions in Eqns. (13T1) and (1331) are fulfilled, corresponding 
to 

01 = - 0s(rin) - 0md(rjn) - (pdipliZin) (42) 

00 = - 0mo(ro«t) (43) 

where rj„ and Tout refer to well inside and outside of the metal surface, respectively. The 
dipole corrections of the K-S potential and the Hellman-Feynman forces in Eqns. (1181) 
and (I22II are now simply obtained by replacing the electrostatic potential 0m (r) + 05(1") 
by the dipole corrected potential 0m (i") + 0s (r) + (pdipo{z) + (pdipi{z) + 0o + 0i. 

Before closing this section, we present the expressions for the double counting terms 
used to evaluate the total energy. In general, the kinetic energy is obtained from the 
one-electron sum and the double counting term as. 



(44) 



Note that adding a constant to the K-S potential, f(r), does not change the kinetic 
energy and 0o + 0i does not give a contribution in this case. The electrostatic part of 
the double counting term is given by. 



Eel = - / Pes(r)[0m(r) + 0s (r) + (l)dipoiz) + (l)dipiiz)]dr (45) 
Jv 

where Pes(i') = —ens{r) is the charge density from the electrons. Adding this term to 
the dipole-corrected, electrostatic potential energy in Eqn. fHTl) . one obtains, 

^ejf,dip ^ ^DC ^ 

\ J A(r)0i(r)rfr- i J Pes(r)0es(r)rfr + (46) 
- pes(r) - Pind{r)]4>dipiiz)dr + (47) 
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^ j Pmd(r)[0i(r) - 0e5(r) + (t)dipi{z)]dY + (48) 
Pi(r)[0mo(r) + (pdipo{z)]dY + (49) 



V 

(00 + \(t>i)Qs. (50) 

where 0i(r) is the ionic potential. Note tliat in the frozen core approximation Pes(r) is 
replaced by the valence charge density and Pi(r) includes the frozen core charge density. 

2.3. Classical perfect conductor model for the metal surface response 

Now we turn to a specific model for the unperturbed charge density Pmo(r) of the 
metal surface and its density response Pmd(r) to an external electrostatic field. In the 
first application of the proposed scheme, we will use the classical perfect conductor (PC) 
model. At the end of this Section, we will discuss how to go beyond this simplistic model 
for the metal surface by using knowledge obtained from previous DFT calculations of 
the semi-infinite jellium model of a metal surface, in which the charge of the nuclei is 
smeared out to a positive homogeneous background. 

The classical PC model is based on a semi-infinite jellium model, here assumed 
to occupy the half-space z < Zp, but makes the bold assumption that the screening 
length is zero, that is, the metallic electrons are able to screen out perfectly any spatial 
variation of the electric field inside the metal. This assumption implies that 

Pmo(r) = (51) 

and also that Pmd(r) will be localised at the jelhum edge at z = Zp, 

Pmd(r) = aind(R)5{z - Zp). (52) 

The surface charge distribution (Tj„d(R) can now be determined by imposing the 
condition of Eqn. (|32l) that the electrostatic potential inside the metal should be zero, 

0(r) = 0, z<Zp, (53) 

and the laterally averaged part of the induced charge density a{z) has to be equal to 

a{z) = a,UG = 0) = (54) 

where Qg is the total charge of S. The laterally varying part of the surface charge 
distribution cr'(r), as obtained from the 2D plane wave coefficients with non-zero 
reciprocal lattice vectors (G ^ 0), can be determined from the condition in Eqn.f l53p . 
Since 0s(r) satisfies the Laplace equation in the metal, the z dependencies of their plane 
wave coefficients for G 7^ are given by, 

(j)s{z, G) = (f)s{zp, G) exp[G'(2; - Zp)] , z < Zp . (55) 

Using the plane wave representation of the Coulomb kernel in Eqn. fl25|) . the 
corresponding coefficients of the induced electrostatic potential from the PC are given 
in this region by, 

4>ind{z, G) = — exp[G{z - Zp)]aind{G) , z < Zp . (56) 
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Using fl55l) the condition in Eqn. fl53p will now be obeyed if and only if, 

G 

CTindiG) = - — (f)s{Zp,G) , (57) 

Ztc 

which in a real space corresponds to the following condition on the laterally varying 
parts of 0s (r) and o"(R), 

^UR) = -^^(R,^p)- (58) 

Note that the derivations of Eqns. ( 1571) and (l58ll are based on the assumption that the 
electron densities of 5* is not overlapping. In the case of such an overlap, the result in 
Eqn.f l58p in contrast to the result in Eqn.( l571) does not obey the symmetry condition for 
the response kernel in Eqn.®. Since in practise a small overlap cannot be avoided we 
will use the result in Eqn.f l57p rather than the result in Eqn. fl58|) . 

The induced electrostatic potential (pindij) can now be obtained outside the PC 
from the induced surface charge distribution using Eqns. flS^ and fl57j) but some care is 
needed to handle the boundary condition for the laterally averaged part of the induced 
potential. The laterally averaged potential 0s (-z) from a charge distribution of S that is 
well-separated from the PC is given by, 

4(^) = ^^^(2; -Zs), z <Zp , (59) 

where Zg is the centroid of the charge distribution of S defined as, 

Ps{z)zdz. (60) 



Qs 

Since according to the condition in Eqn.( 1S^ . (t>ind{zp) = —(f)s{zp), the resulting laterally 
averaged induced electrostatic potential from Eqn.( l54|) in the region outside the PC is 
given by, 

(j^indiz) = ^!I^(^ + Zs- 2Zp) , Z> Zp . (61) 

The plane wave coefficients from the laterally varying part of the induced 
electrostatic potential in the region outside the perfect conductor is now obtained 
directly from Eqns. fl27|) and flFTI) as 

(pindiz, G) = -(t)s{zp, G) exp[-G{z - Zp)] , z > Zp . (62) 

Note that according to Eqns. (I6T]) and ( |62|) . the induced electrostatic potential in real 
space is given outside the PC by the classical image potential corresponding to the 
electrostatic potential from the mirror image of the charge distribution of S in the plane 
z = Zp as given by, 

0i„rf(R, z) = -0s (R, 2zp- z) , z> Zp . (63) 
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2.3.1. Interaction with a point charge In evaluating the proposed scheme based on a 
perfect conductor (PC) model, it is useful to have the result for the interaction energy 
Eint{zo) between the PC and an external point charge Ps(r) = QsS{r — tq), located at 
a position tq in the supercell. The interaction energy is obtained from Eqns.( !TT]) and 
( |4T1) and can be decomposed into one contribution arising from the lateral averaged 
potential, Eint{zo), and one contribution arising from the laterally varying part of the 
potential, E'^^^{zq), 

Eint = ^[0ind(^o) + (Pdiplizo) + 0i] (64) 

EL = Y'^Uzo) (65) 
where, according to Eqns.( l38l) and ( H2l) . the dipole corrections terms are determined by, 

(pdipliz) = ^4—1^0 - Zp)[j- - -\ (66) 

01 = -[(pind{Zp) + (psiZp) + (pdipliZp)]. (67) 

Note that in the absence of perpendicular periodic boundary conditions 4>ind{zo) is given 
by Eqn.(l6T]) and (pdipi{zo) + 0i = 0. The laterally averaged electrostatic potentials are 
given in the supercell by, 

<Ps{z) = ^-[-\z - zo\ + + -] (68) 

(j)ind{z) = —[-\Z-Zp\ + 1^+q\- (69) 

Inserting these results into Eqns.(l67I). one obtains that 

(J^indizo) + (l)dtpliZo) + 01 = ^^^^[{zo - Zp) - ^] , (70) 

and the laterally averaged interaction energy is finally given by, 

EUzo) = ^-^[{zo-Zp)~^] . (71) 

Note that this result differ from the result obtained in the absence of perpendicular 
periodic boundary conditions, 

Eintizo) = ^^^(2:0 - Zp) . (72) 
by the extra term 

^Uzo) = — ^ . (73) 

Furthermore, it is noteworthy that this term diverges when L — t- 00 for fixed A. The 
interaction energy Eint in Eqn.( 172|) is nothing else than the electrostatic energy of a 
parallel plate capacitor where the two plates, each with an area A, are separated by 
the distance zq — Zp and having charges —Qs and Qs- This repulsive interaction energy 
vanishes in the limit A — )■ 00. 
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In the case of the laterally varying part, 4>[ndi^) obtained from ps{z,G) = 
^S{z - zo) and Eqns.(|271) and ([62]) as, 

27tQs ^exp[-G(z + zo-2zp)] r-n fT} m 
0md(r) = 2^ ^ exp [iG.{R - Ro)] 

when the induced potential is well-localised within the supercell corresponding to 
2tt{L — zq)/L\\ ^ 1 and 2t^zq/L\\ ^ 1, with L\\ = L^^y the length of the supercell 
in the x and y directions, parallel to the PC plane. The resulting laterally varying part 
of the interaction energy, E'^^^, in Eqn.( l65|) is then given directly by, 

(.\- 27rg^ ^ exp [-2^(2:0 - Zp)] 
^inA^o) 'a~ ^ 2C ■ ^ ' 

This part of the interaction energy is always attractive and crosses over into the classical 
image potential, 

EU^o) = , , (75) 

A{z - Zo) 

of the point charge in the limit A 00. 

The Hellman-Feynman force Fint{zQ) on the point charge is determined by the dipole 
corrected electric field at the point charge as, 

^int{Zo) = -QsV[0iw(2;o) + (p'indi^o) + (f)dipi{zo) + (^l] (76) 

Using the expressions for (t)ind{zo) and the plane wave coefficients 4>indi'^) 
Eqns.(l6T]). ( l62l) and f l66i) . one obtains directly, 

FL(^o) = ( 1 + E ^M-m^o - z,)]] (77) 

Note that this result for the force is consistent with the result FI^^{zq) = 
—'VEint{zo)z obtained from the interaction energy Eint{zo) = Eint{zo) + E^^^^zq) in 
Eqns.([7I]) and (|74|). 

2.3.2. Beyond the perfect conductor model The first step in going beyond the prefect 
conductor approximation for the metal surface is to account for the non-zero screening 
length of the conduction electrons. Important information and concepts about the 
behaviour of this response have been drawn in the pioneering DFT studies of the semi- 
infinite jellium model based on the LDA by Lang and Kohn[T9]. They showed that there 
is a spill-out of electrons into the vacuum region that creates an extended dipole layer 
with a surface charge density pmo(^)- The dipole moment of this distribution was shown 
to determine the surface contribution to the work function. This charge contribution 
is not accounted for in the PC model. From their calculations of the response of the 
semi-infinite jellium to an external homogeneous electric field, they showed that the 
classical image plane is located at the centroid Zim of the induced density and not at the 
jellium edge as in the perfect conductor model. However, this effect is easily accounted 
for in the PC model by choosing Zp = Zim- 
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The static response of the semi-infinite jeUium to a laterally varying external 
potential can be characterised by a wave- vector dependent reflection coefficient 0(^)121]. 
The 2D plane wave coefficient of the induced electrostatic potential outside the induced 
metal density is given by, 

G) = - g{G)<f)e.t{G) exp[-Gz] (78) 

where (pextiG) is the 2D plane wave coefficient of the external electrostatic potential. 
For example, the interaction energy of a single point charge at z = zq with the metal 
surface becomes, 

E^nt{zo) = -Y J dGexp{-2Gzo)giG) . (79) 



Since the PC model assumes perfect screening for all parallel wave vectors, g{G) = 
exp{2gzp), according to Eqn. (l62|) where Zp gives the position of the PC plane and 
Eqn.(l79ll reduces to the classical image potential in Eqn. (|75l) . Calculations by 
the stabilised jellium model have shown that PC result for g{G) is an excellent 
approximation for G up to G^ =0.8 in the range of 2-4 ao for the electron gas 
density parameter rs[22]. According to Eqn.(179l). this suggests the PC model is still a 
good approximation for the image potential down to distances of about l/2Gc = 0.6 A. 



3. Computational details and implementation 

The proposed scheme that is based on a perfect conductor model has been implemented 
in the plane wave code VASP|23]. The key quantity to compute is the induced 
electrostatic potential 0j„(i(r), which determines how the total energy, K-S potential 
and Hellmann-Feynman forces change in the presence of the PC. The first step is to 
generate aind from the total charge of S, Eqn.f lS^ . and the Fourier components of 
(j^„^(R) from Eqn.( l57|) . The induced density of charge Pindij'^) = o'ind(R)S{z — Zp) was 
then represented at a plane of grid points corresponding to the position of the PC plane, 
from which Pind(g) was generated by the standard routine in VASP based on Eqn.( l35|) . 
The surface dipole moment m that determines (pdip through Eqn.f l38p was obtained from 
Eqn.f l37p . The dipole correction term in Eqn.f l47p is computed from the standard dipole 
correction subroutine by pmd to pes- Finally, the computation of energy term of Eqn.(l48ll 
is carried out in reciprocal space. 

As a simple illustration and test of the proposed scheme, we have considered the 
Na^ ion outside the PC. The electron-ion interaction was described by the projector 
augmented wave method [23]. The six electron in the p semi-core states were treated as 
valence electrons. The electronic exchange and correlation effects were treated within 
the PBE version [25] of the generalised gradient approximation. The plane wave cut-off 
energy was set to the standard value of 400 eV and the Brillouin zone was sampled by 
a 3x3x1 A;-points. 
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Figure 1. Induced surface charge density at the perfect conductor plane for Na"*" 
ion placed at 1.5 A (Upper panel) and 2.0 A (Lower panel) away from the perfect 
conductor. The size of the surface unit cell isloAxlOA and the height of the 
supercell is 20 A. 

4. Results 

The behaviour of the induced surface electron density, (Ti„rf(R), at two different distances 
of the Na+ ion from the perfect conductor plane is shown in Fig. [H for a supercell with 
transversal area A = 100 and a height L of 20 A. It is evident that the induced 
electron density becomes laterally more extended as the ion separates from the surface 
but the net total charge remains constant and is equal to — e. In other words, the 
contributions from plane- wave coefficients cr(G ^ 0) decrease with the distance and 
(jj„rf(R) becomes practically uniform when this distance is sufficiently large. 

The characteristic behaviour of the laterally averaged, dipole-corrected electrostatic 
potential, 4>{z) = 4>s{z) + (pmdiz) + 4>dip{z) + 0i, is illustrated in Fig. [2] as a function of 
the z coordinate when the Na+ ion is placed 2 A outside the perfect conductor plane 
located at z = 7 A. Again, the supercell is the same as in Fig. [H The potential is 
constant and equal to zero inside the perfect conductor as dictated by the condition in 
Eqn.( l53|) . The discontinuous change in slope at the perfect conductor plane (located at 
z = 7A) is induced by the laterally averaged surface charge density a. In the region 
between the perfect conductor plane and the ion (placed ai z = 9A), the slope is 
constant corresponding to a constant electric field. The spatial extension of the charge 
distribution of the ion is reflected by the deviation from this linear behaviour. The 
calculated value of the slope in the potential in this region is 1.81 eV/A, which is 
precisely the value of Ana/eo = 1.81 eV/A for the z-component of the electrical fleld 
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Figure 2. Laterally averaged electrostatic potential, (f>{z) — (f>s{z) + 4>m{z) + (f>dip{z) + 
(/)!, as a function of the z coordinate for a Na+ ion placed 2 A away from the perfect 
conductor plane located at z = 7 A. The dipole layer is located at z = 20 A. Same 
supercell as in Fig. [1] 

inside a capacitor with plates having opposite surface charge densities of a = e/A = 0.01 
eA~^, since A = 100 A^. At the dipole layer located at z = 20 A the potential makes 
a jump, so that the potential is periodic across the supercell boundaries at 2; = and 
20 A. Note that the representation of the surface charge on a single plane of grid points 
give rise to a well-behaved electrostatic potential in contrast to the dipole layer, which 
has to be smeared out in the z direction in order to damp out Gibbs oscillations. 

To facilitate the discussion of the behaviour of the calculated interaction energy, 
Eint, between the Na+ ion and the PC as defined in Eqn.( !TT|) . we have decomposed 
this energy into the contributions Eint and from the laterally averaged, amd, and 
the laterally varying surface charge density, (jj'^^(R), respectively. Note that EsIusq] 
in Eqn. lfTT]) has been obtained for the ion in the supercell with a uniform neutralising 
background. The two contributions Ei^t and i^^'^j are shown in Figs. [3]and|l]for different 
surface areas of the supercell. Since we shall compare DFT results for the Na^ with the 
analytical results for a point charge in Eqns. ( |7T]l and (Clj), we need to avoid a significant 
overlap between the electron density of the ion and the induced electron density at the 
perfect conductor plane. With this, we shall only consider distances of the Na"*" from 
the perfect conductor plane larger than 0.8 A. 

As shown in Fig. |3] (Upper panel), the contribution from the uniform surface 
charge distribution results in a linear variation of the interaction energy with the ion 
distance to the PC surface. This linear behaviour and the decrease in magnitude of 
the associated force (given by the negative gradient) with increasing transversal area 
is in agreement with the corresponding result of Eqn.( !7T|) in the point charge model. 
The relative difference between the latter result and the computed interaction energy is 
smaller than the 2%, as shown in Fig. |3] (Lower panel). These small differences can in 
part be attributed to the polarisation of the semi-core of the Na"'"and at shorter distances 
in part the spatial extension of the charge distribution of the ion. Note that E^nt for 
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Figure 3. (Upper panel) Calculated contribution, Eint, to the interaction energy from 
the laterally averaged part of the induced surface charge density, ct, as a function of 
the distance to the perfect conductor plane for the different sizes of the surface unit 
cell. (Lower panel) Calculated energy difference with respect to the point charge case 
result of Eqn.dni)- The height of the superceh is 20 A. 



the different A do not cross at the perfect conductor plane but at a distance of about 
1.75 A outside the perfect conductor plane due to the extra energy term of Eqn. (173l) . In 
fact this term gives that the crossing is located at L/12, which is in excellent agreement 
with the result in Fig. |3] (Upper panel). 

In contrast to the laterally average interaction energy Ei^t, the contribution to the 
interaction energy from the laterally varying part E'^^^ is always attractive, becomes more 
dominant with increasing surface area and decay rapidly with the distances from the 
perfect conductor, as shown in Fig. H] (Upper panel). This result is in close agreement 
with the result in Eqn.f TMl) . obtained from the corresponding contribution in the point 
charge model, as shown by the difference between these two results in Fig. H] (Lower 
panel). In fact, relative energy differences between DFT and the point charge case 
are smaller than 2%. In the limit of infinite surface area, the interaction tends to the 
classical image interaction. 

We now present the calculated Hellmann-Feynman force along the z direction in 
Fig. [5l The force is always attractive and becomes constant for large distances of the ion 
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Figure 4. (Upper panel) Calculated contribution, i?j'„j, to the interaction energy from 
the laterally varying part of the induced surface charge density, (t'(R), as a function 
of the distance of the ion to the perfect conductor plane for the different sizes of the 
surface unit cell. (Lower panel) Calculated energy difference with respect to the point 
charge model of Eqn. ([74l) . The height of the supercell is 20 A. 

from the perfect conductor plane. As the ion approaches to perfect conductor, the effect 
of the laterally varying components of the induced charge density become increasingly 
more important, thus leading to an increased attraction. Clearly, the magnitude of this 
attraction increases with the transversal area. In the limit of infinite surface area, the 
induced force will tend to the classical force exerted by the image charge. We observe 
that the differences with respect to the point charge case are smaller than 2.5%. 

Finally, as a critical test of the implementation of the perfect conductor model, 
we analyse the consistency between the computed Hellman-Feynman force on the ion 
and the gradient of the total energy. In Fig. El we show the difference between the 
force, Fz, computed along the z direction using Eqn. (l22|) and the numerical derivative 
of the interaction energy, —dEint/dzo, as a function of the Na^ distance to the perfect 
conductor. Ideally, this difference should be zero but, due to numerical errors, small 
deviations are always expected. In fact, errors are essentially smaller than the 0.06% of 
the computed forces, showing that sufficient numerical consistency has been achieved. 
This level of self-consistency is also obtained for smaller distances where we have a more 
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Figure 5. Calculated Hellman-Feynman force on the Na+ ion as a function the 
distance, zq, from the perfect conductor plane for the different sizes of the surface unit 
cell. The height of the super cell is 20 A. Differences with respect to the point charge 
case are lower than 2.5%. 
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Figure 6. Test of the consistency between the Hellman-Feynman force on the ion and 
the force from the gradient of the total energy. AF is the difference between these two 
forces and values are smaller than the 0.06% of the computed forces. Same supercell 
as in Fig. [T] 



substantial overlap of the electron densities. Note that this high degree of consistency 
even in the case of overlapping electron densities is guaranteed by using the density 
response in Eqn.f l57|) that obeys the symmetry condition in Eqn.®. 

The above results suggest that we have derived and implemented a DFT method for 
a charged system in front of a perfect conductor. Small differences with respect to the 
the point charge case could be attributed to the differences in the spatial extension of the 
densities of charge and that the point charge model does not include any polarisation, 
which will be inevitably present when placing atoms or molecules under the electrical 
field generated by the induced potential. 

We believe that the use of this new DFT methodology becomes particularly 
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convenient since it allows to the possibility of computing DFT problems of different 
charge states in a controllable way, that is, by defining (at will) the amount of charge 
transfer between the system and the perfect conducting plane. However, we have shown 
that the use of the perfect conductor approximation only induces attractive interactions, 
which would make the system to move closer and closer to the perfect conductor plane. 
This represent a limitation if one aims to make use of this DFT methodology to simulate 
realistic problems involving metallic surfaces. In fact, the metallic surface will exert 
repulsive forces for sufficiently close distances, thus avoiding the system to collapse with 
the surface [26]. This lack of repulsion is a direct consequence of having considered 
only the electrostatic interactions in the approximate energy functional of Eqn.(Il]). In 
a following publication [27J, we propose a new procedure to surmount this limitation. 

5. Concluding Remarks 

To the final purpose of developing a simplified density functional theory (DFT) method 
for treating charged atoms and molecules on an ultrathin-insulating film supported by 
a metal substrate, we have presented a new approximate DFT methodology for the 
calculation of the total energy, ionic (Hellman-Feynman) forces and electronic structure 
of a charged system placed in front of a metal surface. In this new methodology, 
the electron densities of the metal surface and the charged system are assumed to 
be non-overlapping and there is only an electrostatic interaction between these two 
fragments. Whereas the charged system is treated fully using DFT, the metal surface 
is approximated within linear response, corresponding to an expansion of the energy 
of the metal surface to second order in the electrostatic potential from the charged 
system. In particular, we have carried out a careful analysis of the effect of periodic 
boundary conditions and derived appropriate dipole corrections for the total energy 
and the ionic forces. The proposed method have two main advantages. First, the 
metal electron states are not treated explicitly but only implicitly through their density 
response to an external electrostatic field, which can be captured in a simple model. 
Here, we have explicitly considered the perfect conductor approximation for the metal 
response although the method is not strictly limited to this approximation. In the 
perfect conductor model the screening length by the conduction electron is assumed to 
be zero and the screening charge resides on a plane. Second, different charge states 
of the charged system can be handled directly. Based on the simple perfect conductor 
approximation for the metal response, we have implemented this method in the density 
functional theory VASP code. 

A simple illustration and test of this method is provided by the case of the Na"*" 
ion outside a perfect conductor. The six electrons in the p semi-core of the ion were 
included in the calculations. The success of our implementation is demonstrated by the 
excellent agreement between the calculated Hellman-Feynman force and the gradient of 
the interaction energy, even in the region close to the perfect conductor plane where 
the electron density of the ion overlaps with the induced charge density at the perfect 
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conductor plane. Furthermore, the small polarisability of the ion makes the calculated 
interaction energy to be in close agreement to the analytical result for a point charge 
outside the perfect conductor within a supercell of finite size. 

Finally, to fulfil our overall aim, we need to include the missing interactions in our 
simplified scheme arising from over-lapping densities and van der Waals interactions 
between the ultrathin-insulating film with charged adsorbates and the metal substrate. 
We plan to do this by developing a simple parametrised force field between the metal 
substrate and the atoms in the adjacent layer of the film, where the parameters are 
obtained by fitting the resulting interactions to DFT calculations of the film (without 
adsorbates) over the metal substrate. The presentation of this scheme and associated 
results of this model are deferred to a separate publication. 
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